/*
 * OdeSolver.h
 *
 *  Created on: 28 Jul 2011
 *      Author: Allan
 */

#ifndef ODESOLVER_H_
#define ODESOLVER_H_

// Eigen includes
#include <Eigen/Dense>
using namespace Eigen;

// GSL includes
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

class OdeProblem
{
public:
	virtual void Function(double t, const VectorXd& y, VectorXd& f) = 0;

	virtual void Jacobian(double t, const VectorXd& y, MatrixXd& dfdy, VectorXd& dfdt) {};
};

class OdeSolver
{
public:
	enum Algorithm
	{
		/// The list of supported algorithms (for more information, check: http://www.gnu.org/software/gsl/manual/html_node/Stepping-Functions.html)
		RungeKutta2,                 // GSL name: (gsl_odeiv2_step_rk2)
		RungeKutta4,                 // GSL name: (gsl_odeiv2_step_rk4)
		RungeKuttaFehlberg45,        // GSL name: (gsl_odeiv2_step_rkf45)
		RungeKuttaCashKarp45,        // GSL name: (gsl_odeiv2_step_rkck)
		RungeKuttaPrinceDormand8,    // GSL name: (gsl_odeiv2_step_rk8pd)
		RungeKuttaImplicitGaussian1, // GSL name: (gsl_odeiv2_step_rk1imp)
		RungeKuttaImplicitGaussian2, // GSL name: (gsl_odeiv2_step_rk2imp)
		RungeKuttaImplicitGaussian4, // GSL name: (gsl_odeiv2_step_rk4imp)
		BulirschStoerImplicit,       // GSL name: (gsl_odeiv2_step_bsimp)
		MultistepAdams,              // GSL name: (gsl_odeiv2_step_msadams)
		MultistepBDF                 // GSL name: (gsl_odeiv2_step_msbdf)
	};
	
public:
	OdeSolver();
	
	OdeSolver(unsigned dimension);
	
	~OdeSolver();
	
	void
	Initialiaze(unsigned dimension);
	
	void
	SetAlgorithm(Algorithm algorithm);
	
	void
	SetControlStandard(double epsAbs, double epsRel, double yFactor, double dydtFactor);
	
	void
	SetControlSolution(double epsAbs, double epsRel);
	
	void
	SetControlSolutionDerivative(double epsAbs, double epsRel);
	
	bool
	Iterate(OdeProblem& problem, double& t, double tEnd, double& dt, VectorXd& y);
	
	bool
	Solve(OdeProblem& problem, double& t, double tEnd, double& dt, VectorXd& y);

public:
	struct Params
	{
		OdeProblem* problem;
		VectorXd y;
		VectorXd f;
		VectorXd dfdt;
		MatrixXd dfdy;
	};
	
private:
	unsigned dimension;
	
	Params params;
	
	const gsl_odeiv2_step_type* type;
	
	gsl_odeiv2_step* step; 
	
	gsl_odeiv2_control* control;
	
	gsl_odeiv2_evolve* evolve;
	
	double* yGSL;
};

#endif /* ODESOLVER_H_ */

